% code to compute MT trfn from raw ts files !!!
% written for svs
% basically to yeild tfs at all FC's
% clears unneceesary vriables after
% invoking readrts reads the header
% compute cal. table according to gain settings etc;
% computes cross and auto spectra & average it 
% MT app rest & phse are now computed
% - Bug report: fixed the conjugate problem
% latest date 6.10.5

readrts;

Notchi = NotchF(2:5);		% NOTE ! there may be problems 
IsNotch = Filt(2:5);		%        if notch do not have right dimension..
hpb1 = Filt(1);			% NOTE ! high pass filter is always disabled
prob = [Exd Eyd Eyd Eyd Eyd];	% 	 in calib05.m
gain = [Exa Eya Hxa Hya Hza];
clear Bden Record3 Bdend Bddm Count Record4 hpf Bdst Record5;
clear Bdstart MaxRecord ans nt Filt NBlocks array nlines Events ;
clear NotchF dummy Record1 factor Channel;
BandFrq = [8192 256 8 1/4 1/128 1/4096]; % Here the upper and lower limit of frequencies of each band is given
RadSize = [0.125 0.25 0.5 0.75];
choiseR = 1;
block = 1024;
TLFreq = logspace(log10(BandFrq(nband)),log10(BandFrq(nband+1)),20);
TLFreq([1 end]) = [];
notch = Notchi(find(IsNotch>0)); 
PAI = 3.1415926535897;

ProcDef = setdefault(nband);

[TLine,TLRad,TLFreq,TBand] =  SetTLine(block,RadSize(choiseR));
ProcDef.TLFreq1 = TLFreq;
for i = 1:nchannel,
[cplx,frq1] = calib05(i,nband,block,prob(i),notch,hpb1,gain(i));
CalData(i,:) = cplx;
end;

%------------native declarations---------------

WI = hanning(block);
%len = length(matrix)/block;
len=50;
start = 0;
st = start*block + 1;
en = (start+1)*block;
nchannel = 5;


nfrq = length(frq1);
HxHx = zeros([nfrq,1]);

HyHy = HxHx;
HxHy = HxHx;
HyHx = HxHx;
ExHx = HxHx;
EyHx = HxHx;
ExHy = HxHx;
EyHy = HxHx;


for i = (start+1):len,

%fft & calib
for u = 1:5,
	K = CorTre(matrix(st:en,u),0);
	Channel1(:,u) = fft(K.*WI)*2/block;
	Channel(:,u) = Channel1(2:(block/2)+1,u).*conj(CalData(u,1:block/2)');
end; 	
   
   Hx   = Hx + Channel(:,3);
   Hy   = Hy + Channel(:,4);
   Hz   = Hz + Channel(:,5);
   Ex   = Ex + Channel(:,1);
   Ey   = Ey + Channel(:,2);
%------------------------------------------------

%clear Channel;

st = en + 1;
en = (i+1)*block;
end;


Hy =Hy/len;
Hx =Hx/len;
Ex = Ex/len;
Ey = Ey/len;
Hz = Hz/len;


